knitr::opts_chunk$set(echo = TRUE)
library("ggplot2")
library("R.matlab")
library("reshape2") 
library("pracma") 
library("plot.matrix")
library("MASS")
library("ggcorrplot")
library("cowplot")
library("corrplot")
library("Rfast")
library("factoextra")
library("gridExtra")

# Global variables
N <- 240
V <- 441
x1 <- 21
x2 <- 21
nsrcs <- 6

Question 1.1

# function that generate TCs
generateTC <- function(ones_onsets, DO, timeseries) 
{
    stim <- numeric(timeseries)
    for(i in 1:length(ones_onsets))
    {
        stim[(ones_onsets[i]+1):(DO+ones_onsets[i])] <- c(replicate(DO,1))
    }
    return(stim)
}
AV <- c(0, 20, 0, 0, 0, 0)
IV <- c(30, 45, 60, 40, 40, 40)
DO <- c(15, 20, 25, 15, 20, 25)
TC <- matrix(0, nrow = 240, ncol = 6)
# Generate TCs
for (i in 1:6) {
  TC[,i] <- generateTC(seq(AV[i], N, by = IV[i]), DO[i], N)[1:N]
}

TC <- data.frame(TC)
par(mfrow=c(3,2))

# Standardize TC and plot TCs
for (i in 1:6)
{
  TC[,i] <- (TC[,i] - mean(TC[,i])) / sd(TC[,i]);
  plot(TC[, i], typ="l", main=paste("TC", i), ylab=paste("Temporal Source", i))
}

Question 1.2

corrplot(cor(TC), method = 'color', order = 'alphabet', title="\nCorrelation Matrix", type = "lower", tl.srt = 30)

Question 1.3

# Construct tmpSM
tmpSM <- array(0, c(6, 21, 21))
tmpSM[1, 2:6, 2:6] = 1
tmpSM[2, 2:6, 15:19] = 1
tmpSM[3, 8:13, 2:6] = 1
tmpSM[4, 8:13, 15:19] = 1
tmpSM[5, 15:19, 2:6] = 1
tmpSM[6, 15:19, 15:19] = 1

# Plot SMs in six subplots
par(mfrow=c(2, 3))
for (i in 1:6)
{
  plot(tmpSM[i,,], main=paste("tmpSM",i), xlab=NA, ylab=NA, border=NA)
}

# plot 6 vectored SMs to show independence
SM <- matrix(c(c(tmpSM[1,,]), c(tmpSM[2,,]), c(tmpSM[3,,]), c(tmpSM[4,,]), c(tmpSM[5,,]), c(tmpSM[6,,])), 6, 21*21, byrow=TRUE)

corrplot(cor(t(SM)), main="\nSM Correlation Heatmap", method="color")

Question 1.4

# Construct Gamma t and Gamma s
set.seed(10)
Gamma_t = matrix(rnorm(240*6, mean=0, sd=0.25**0.5), 240, 6)
Gamma_s = matrix(rnorm(6*441, mean=0, sd=0.015**0.5), 6, 441)

# Plot correlation
par(mfrow=c(1, 2))
corrplot(cor(Gamma_t), main="\n\n\n\nTemporal Noise", 
         xlab=NA, ylab=NA, method="color")
corrplot(cor(t(Gamma_s)), main="\n\n\n\nSpatial Noise", 
         xlab=NA, ylab=NA, method="color")

df = data.frame(c(Gamma_t))
# Plot histogram of both noise source
gamma_t_plot <- ggplot(df, aes(c(c.Gamma_t.)))+
  geom_histogram(aes(y=..density..),
  binwidth=0.15, color="darkblue", fill="lightblue") +
  labs(title="Distribution of Temporal Noise",x="Value", y = "Density") +
  stat_function(fun = dnorm, args = list(mean = 0, sd = (0.25)**0.5)) +
  geom_vline(xintercept = (1.96)*(0.25)**0.5, color="red") + 
  geom_vline(xintercept = -(1.96)*(0.25)**0.5, color="red")

df = data.frame(c(Gamma_s))
gamma_s_plot <-ggplot(df, aes(c(c.Gamma_s.)))+
  geom_histogram(aes(y=..density..),
  binwidth=0.05, color="darkblue", fill="lightblue") +
  labs(title="Distribution of Spatial Noise",x="Value", y = "Density") +
  stat_function(fun = dnorm, args = list(mean = 0, sd = (0.015)**0.5)) +
  geom_vline(xintercept = (1.96)*(0.015)**0.5, color="red") + 
  geom_vline(xintercept = -(1.96)*(0.015)**0.5, color="red")

plot_grid(gamma_t_plot, gamma_s_plot)

Gamma_ts <- Gamma_t%*%Gamma_s
# Check correlation of the product of Gamma t and Gamma s
corr <- cor(Gamma_ts)

# Randomly sample 10 variables to visualise
samples <- sample (c(1:441), size=10, replace=F)
corrplot(corr[samples, samples], method="color", main="\nCorrelation between 441 variables (Sampled)", type = "lower", tl.srt = 30)

Question 1.5

# Generate a synthetic dataset X
TC <- as.matrix(TC, 6, 240)
X <- (TC + Gamma_t) %*% (SM + Gamma_s)

# Check compatibility
c(dim(TC), dim(Gamma_s))
[1] 240   6   6 441
c(dim(Gamma_t), dim(SM))
[1] 240   6   6 441
dim(X)
[1] 240 441
# Plot at least 100 randomly selected time-series from X
samples <- data.frame(n=1:240, X[, sample.int(240, 100)])
samples <- melt(samples, id.vars = "n")
ggplot(samples, aes(x=n, y=value, col=variable)) + geom_line() + 
  labs(title="100 Ramdonly Selected Time-series", x="", y="Value") +
  theme(plot.margin = margin(1.5,.5,1.5,.5, "cm"),
        legend.key.height= unit(0.1, 'cm'),
        legend.key.width= unit(0.1, 'cm'))

# plot variance of all 441 variables
variances <- colVars(X)
variances <- data.frame(n=1:441, variances)
ggplot(variances, aes(x=n, y=variances)) +geom_point(color="steelblue") + 
  labs(title="Variance of 441 variables", x="Variables", y="Variances")

# Standardize X
X = scale(X)

Question 2.1

D <- as.matrix(TC)

# Estimate A using least square estimation
A_lsr <- abs(solve(t(D)%*%D)%*%t(D)%*%X)

# Calculate D_lsr
D_lsr <- X%*%t(A_lsr)
pal <- colorRampPalette(c("blue", "yellow"))
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), 6, 2, ncol = 2), 
       widths=c(2.5,5))
par(mar = c(2,5,1.5,1.5))
for (i in 1:6) 
{
  A_lsr_p <- A_lsr[i,]
  D_lsr_p <- D_lsr[,i]
  dim(A_lsr_p) <- c(21, 21)
  plot(A_lsr_p, border=F, col = pal(40), xlab=NA, ylab=NA, 
       main=paste("Retrieved Spatial Map",i))
  
  plot(D_lsr_p, type = "l", xlab=NA, ylab=NA, 
       main=paste("Retrieved Time Course",i))
}

df <- data.frame(D_lsr[,3], X[,30])
colnames(df) <- c("D_lsr_3", "X_30")
plot_1 <- ggplot(df, aes(y=D_lsr_3, x=X_30)) + geom_point(color = "steelblue") +
  labs(title=expression(paste(D[lsr], "  column 3 vs X column 30")),
       x=expression(paste(D[lsr], " column 3")), y = "X column 30") + geom_smooth()

df <- data.frame(D_lsr[,4], X[,30])
colnames(df) <- c("D_lsr_4", "X_30")
plot_2 <- ggplot(df, aes(y=D_lsr_4, x=X_30)) + geom_point(color = "steelblue") +
  labs(title=expression(paste(D[lsr], "  column 4 vs X column 30")),
       x=expression(paste(D[lsr], " column 4")), y = "X column 30") + geom_smooth()

suppressMessages(plot_grid(plot_1, plot_2))

Question 2.2

lambda <- 0.5
I <- matrix(0, 6, 6)
diag(I) <- 1
# Compute scalar value
penalty_term <- lambda * V
# Estimate A_rr and D_rr
A_rr <- abs(solve(t(D)%*%D + penalty_term*I)%*%t(D)%*%X)
D_rr <- X%*%t(A_rr)

C_tlsr = 0
C_trr = 0
for (i in 1:6)
{
  C_tlsr[i] = cor(TC[,i], D_lsr[,i])
  C_trr[i] = cor(TC[,i], D_rr[,i])
}

# Compare the sum of these two correlation vectors
print(paste("Sum of C_tlsr:", sum(C_tlsr)))
[1] "Sum of C_tlsr: 5.09748114621307"
print(paste("Sum of C_trr:", sum(C_trr)))
[1] "Sum of C_trr: 5.16285782879284"
# Set lambda value as 1000 for plot visualisation
lambda <- 1000
# Compute scalar value
penalty_term <- lambda * V
# Estimate A_rr and D_rr
A_rr_1000 <- abs(solve(t(D)%*%D + penalty_term*I)%*%t(D)%*%X)

# plot and compare the shrinkage of variables
values <- data.frame(n=1:441, val=A_lsr[1,])
A_lsr_plot <- ggplot(values, aes(x=n, y=val)) +
  geom_point(color = "steelblue") + 
  labs(title=expression(A[lsr]), y="Values", x="Index")

values <- data.frame(n=1:441, val=A_rr_1000[1,])
A_rr_plot <- ggplot(values, aes(x=n, y=val)) +
  geom_point(color = "steelblue") + 
  labs(title=expression(paste(A[rr], "   ", lambda, " = 1000")), y="Values",
       x="Index")

suppressMessages(plot_grid(A_lsr_plot, A_rr_plot))

Question 2.3

# Select ρ between 0 and 1 with an interval of 0.05
rho_values <- seq(from=0, to=1, by=0.05)
MSE <- rep(0, 21)
for (j in 1:21)
{
  rho <- rho_values[j]
  step <- 1/(norm(TC %*% t(TC)) * 1.1)
  thr <- rho*N*step
  Ao <- matrix(0, nsrcs, 1)
  A <- matrix(0, nsrcs, 1)
  Alr <- matrix(0, nsrcs, x1*x2)
  MSE_temp <- rep(0, 10)
  
  # Generate new standardized X in terms of Gamma t and Gamma s
  for (a in 1:10) {
    set.seed(a)
    Gamma_t_temp = matrix(rnorm(240*6, mean=0, sd=0.25**0.5), 240, 6)
    Gamma_s_temp = matrix(rnorm(6*441, mean=0, sd=0.015**0.5), 6, 441)
    X_temp <- (TC + Gamma_t_temp) %*% (SM + Gamma_s_temp)
    X_temp <- scale(X_temp)
    
    for (k in 1:(x1*x2)) {
      A <- Ao+step*(t(TC) %*% (X_temp[,k]-(TC%*%Ao)))
      A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
      
      for (i in 1:10) {
        Ao <- A
        A <- Ao+step * (t(TC)%*%(X_temp[,k]-(TC%*%Ao)))
        A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
        }
      Alr[,k] <- A
    }
    Alr <- Alr
    Dlr <- X_temp%*%t(Alr)
    MSE_temp[a] <- norm((X_temp - Dlr%*%Alr), type='F')^2/(N*V)
  }
  
  # Compute the average MSE of the current ρ
  MSE[j] <- mean(MSE_temp)
}
# plot average of MSE over these 10 realizations against each value of ρ

df <- data.frame(rho_values, MSE)
plot_1 <- ggplot(df, aes(x=rho_values, y=MSE)) + geom_point(color="steelblue") +
  labs(title="Average MSE for each ρ", x="ρ values", y="Average MSE")

df <- df[10:21,]
plot_2 <- ggplot(df, aes(x=rho_values, y=MSE)) + geom_point(color="steelblue") +
  labs(title="Average MSE for each ρ (Zoomed In)", x="ρ values", y="Average MSE")

suppressMessages(plot_grid(plot_1, plot_2))

Question 2.4

# Estimate A_lr and D_lr using the selected ρ value
rho <- 0.6
step <- 1/(norm(TC %*% t(TC)) * 1.1)
thr <- rho*N*step
Ao <- matrix(0, nsrcs, 1)
A <- matrix(0, nsrcs, 1)
A_lr <- matrix(0, nsrcs, x1*x2)

for (k in 1:(x1*x2)) {
  A <- Ao+step*(t(TC) %*% (X[,k]-(TC%*%Ao)))
  A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  
  for (i in 1:10) {
    Ao <- A
    A <- Ao+step * (t(TC)%*%(X[,k]-(TC%*%Ao)))
    A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  }
  A_lr[,k] <- A
}
A_lr <- abs(A_lr)
D_lr <- X%*%t(A_lr)
# estimate four correlation vectors
C_trr <- 0
C_srr <- 0
C_tlr <- 0
C_slr <- 0

for (i in 1:6)
{
  C_trr[i] = cor(TC[,i], D_rr[,i])
  C_tlr[i] = cor(TC[,i], D_lr[,i])
  C_srr[i] = cor(SM[i,], A_rr[i,])
  C_slr[i] = cor(SM[i,], A_lr[i,])  
}
sum(C_trr)
[1] 5.162858
sum(C_tlr)
[1] 5.294203
sum(C_srr)
[1] 3.388821
sum(C_slr)
[1] 5.043101
layout(matrix(seq(from=1, to=24, by=1), 6, 4, ncol = 4), 
       widths=c(3,5, 3, 5))
par(mar = c(3,5,1.5,2))
for (i in 1:6) {
  A_rr_p <- A_rr[i,]
  A_lr_p <- A_lr[i,]
  dim(A_rr_p) <- c(21, 21)
  dim(A_lr_p) <- c(21, 21)
  plot(A_rr_p, border=F, col = pal(40), xlab=NA, ylab=NA, 
       cex.main=1.6, cex.lab=1.5, cex.axis=1.2, main=paste("Arr",i))
  plot(D_rr[,i], xlab=NA, ylab=NA, main=paste("Drr",i), 
       cex.main=1.6, cex.lab=1.5, cex.axis=1.2, col="blue", type="l")
  plot(A_lr_p, border=F, col = pal(40), xlab=NA, ylab=NA, 
       cex.main=1.6, cex.lab=1.5, cex.axis=1.2, main=paste("Alr",i))
  plot(D_lr[,i], xlab=NA, ylab=NA, main=paste("Dlr",i),
       cex.main=1.6, cex.lab=1.5, cex.axis=1.2, col="blue", type="l")
}

Question 2.5 using SVD


pcr <- svd(D)
Z <- pcr$u
par(mfrow=c(6, 2), mar = c(2,3,1.5,2))

# Plot the regressors in Z and source TCs side by side
for (i in 1:6) {
  plot(TC[,i], type="l", main=paste("TC", i), cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
  plot(Z[,i], type="l", main=paste("Z", i), cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
}

ev <- pcr$d
ev <- data.frame(index=seq(1, 6), ev)
plot(ev, col="blue", main="Eigenvalues of each PC", xlab="PCs", ylab="Eigenvalues",
     cex.main=1.4, cex.lab=1.4, cex.axis=1.4)

# Estimate A_pcr and D_pcr using the ρ = 0.001
matrix <- Z
rho <- 0.001
step <- 1/(norm(Z %*% t(Z)) * 1.1)
thr <- rho*N*step
Ao <- matrix(0, nsrcs, 1)
A <- matrix(0, nsrcs, 1)
A_pcr <- matrix(0, nsrcs, x1*x2)

for (k in 1:(x1*x2)) {
  A <- Ao+step*(t(matrix) %*% (X[,k]-(matrix%*%Ao)))
  A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  
  for (i in 1:10) {
    Ao <- A
    A <- Ao+step * (t(matrix)%*%(X[,k]-(matrix%*%Ao)))
    A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  }
  A_pcr[,k] <- A
}
A_pcr <- abs(A_pcr)
D_pcr <- X%*%t(A_pcr)
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), 6, 2, ncol = 2), 
       widths=c(0.8,2))
par(mar = c(2,5,1.5,1.5))

# Plot the results of DPCR and APCR side by side
for (i in 1:6) 
{
  A_pcr_p <- A_pcr[i,]
  D_pcr_p <- D_pcr[,i]
  dim(A_pcr_p) <- c(21, 21)
  plot(A_pcr_p, border=F, col = pal(40), xlab=NA, ylab=NA, 
       main=paste("A pcr",i))
  
  plot(D_pcr_p, type = "l", xlab=NA, ylab=NA, 
       main=paste("D pcr",i))
}

Question 2.5 using PRCOMP

m <- prcomp(D, center = TRUE, scale.=TRUE)
Z <- m$x
pev <- fviz_eig(m) + 
  labs(title="Percent of explained variances of PCs", x="PCs")
eigenvalues <- fviz_eig(m, choice="eigenvalue", geom="line") + 
  labs(title="Eigenvalue of each PC", x="PCS")
suppressMessages(plot_grid(pev, eigenvalues))

par(mfrow=c(6, 2), mar = c(2,3,1.5,2))

# Plot the regressors in Z and source TCs side by side
for (i in 1:6) {
  plot(TC[,i], type="l", main=paste("TC", i), cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
  plot(Z[,i], type="l", main=paste("Z", i), cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
}

# Estimate A_pcr and D_pcr using the ρ = 0.001
matrix <- Z
rho <- 0.001
step <- 1/(norm(Z %*% t(Z)) * 1.1)
thr <- rho*N*step
Ao <- matrix(0, nsrcs, 1)
A <- matrix(0, nsrcs, 1)
A_pcr <- matrix(0, nsrcs, x1*x2)

for (k in 1:(x1*x2)) {
  A <- Ao+step*(t(matrix) %*% (X[,k]-(matrix%*%Ao)))
  A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  
  for (i in 1:10) {
    Ao <- A
    A <- Ao+step * (t(matrix)%*%(X[,k]-(matrix%*%Ao)))
    A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  }
  A_pcr[,k] <- A
}
A_pcr <- abs(A_pcr)
D_pcr <- X%*%t(A_pcr)
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), 6, 2, ncol = 2), 
       widths=c(0.8,2))
par(mar = c(2,5,1.5,1.5))

# Plot the results of DPCR and APCR side by side
for (i in 1:6) 
{
  A_pcr_p <- A_pcr[i,]
  D_pcr_p <- D_pcr[,i]
  dim(A_pcr_p) <- c(21, 21)
  plot(A_pcr_p, border=F, col = pal(40), xlab=NA, ylab=NA, 
       main=paste("A pcr",i))
  
  plot(D_pcr_p, type = "l", xlab=NA, ylab=NA, 
       main=paste("D pcr",i))
}

df <- data.frame(m$rotation)
df
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkKbGlicmFyeSgiZ2dwbG90MiIpCmxpYnJhcnkoIlIubWF0bGFiIikKbGlicmFyeSgicmVzaGFwZTIiKSAKbGlicmFyeSgicHJhY21hIikgCmxpYnJhcnkoInBsb3QubWF0cml4IikKbGlicmFyeSgiTUFTUyIpCmxpYnJhcnkoImdnY29ycnBsb3QiKQpsaWJyYXJ5KCJjb3dwbG90IikKbGlicmFyeSgiY29ycnBsb3QiKQpsaWJyYXJ5KCJSZmFzdCIpCmxpYnJhcnkoImZhY3RvZXh0cmEiKQoKIyBHbG9iYWwgdmFyaWFibGVzCk4gPC0gMjQwClYgPC0gNDQxCngxIDwtIDIxCngyIDwtIDIxCm5zcmNzIDwtIDYKYGBgCiMgUXVlc3Rpb24gMS4xCmBgYHtyfQojIGZ1bmN0aW9uIHRoYXQgZ2VuZXJhdGUgVENzCmdlbmVyYXRlVEMgPC0gZnVuY3Rpb24ob25lc19vbnNldHMsIERPLCB0aW1lc2VyaWVzKSAKewoJc3RpbSA8LSBudW1lcmljKHRpbWVzZXJpZXMpCglmb3IoaSBpbiAxOmxlbmd0aChvbmVzX29uc2V0cykpCgl7CgkJc3RpbVsob25lc19vbnNldHNbaV0rMSk6KERPK29uZXNfb25zZXRzW2ldKV0gPC0gYyhyZXBsaWNhdGUoRE8sMSkpCiAgIAl9CglyZXR1cm4oc3RpbSkKfQpgYGAKCmBgYHtyfQpBViA8LSBjKDAsIDIwLCAwLCAwLCAwLCAwKQpJViA8LSBjKDMwLCA0NSwgNjAsIDQwLCA0MCwgNDApCkRPIDwtIGMoMTUsIDIwLCAyNSwgMTUsIDIwLCAyNSkKVEMgPC0gbWF0cml4KDAsIG5yb3cgPSAyNDAsIG5jb2wgPSA2KQojIEdlbmVyYXRlIFRDcwpmb3IgKGkgaW4gMTo2KSB7CiAgVENbLGldIDwtIGdlbmVyYXRlVEMoc2VxKEFWW2ldLCBOLCBieSA9IElWW2ldKSwgRE9baV0sIE4pWzE6Tl0KfQoKVEMgPC0gZGF0YS5mcmFtZShUQykKYGBgCgpgYGB7ciwgZWNobz1UUlVFLCBmaWcuaGVpZ2h0PTMuNSwgZmlnLndpZHRoPTMuOH0KcGFyKG1mcm93PWMoMywyKSkKCiMgU3RhbmRhcmRpemUgVEMgYW5kIHBsb3QgVENzCmZvciAoaSBpbiAxOjYpCnsKICBUQ1ssaV0gPC0gKFRDWyxpXSAtIG1lYW4oVENbLGldKSkgLyBzZChUQ1ssaV0pOwogIHBsb3QoVENbLCBpXSwgdHlwPSJsIiwgbWFpbj1wYXN0ZSgiVEMiLCBpKSwgeWxhYj1wYXN0ZSgiVGVtcG9yYWwgU291cmNlIiwgaSkpCn0KYGBgCiMgUXVlc3Rpb24gMS4yCmBgYHtyfQpjb3JycGxvdChjb3IoVEMpLCBtZXRob2QgPSAnY29sb3InLCBvcmRlciA9ICdhbHBoYWJldCcsIHRpdGxlPSJcbkNvcnJlbGF0aW9uIE1hdHJpeCIsIHR5cGUgPSAibG93ZXIiLCB0bC5zcnQgPSAzMCkKYGBgCiMgUXVlc3Rpb24gMS4zCmBgYHtyLCBlY2hvPVRSVUUsIGZpZy5oZWlnaHQ9Mi41LCBmaWcud2lkdGg9NH0KIyBDb25zdHJ1Y3QgdG1wU00KdG1wU00gPC0gYXJyYXkoMCwgYyg2LCAyMSwgMjEpKQp0bXBTTVsxLCAyOjYsIDI6Nl0gPSAxCnRtcFNNWzIsIDI6NiwgMTU6MTldID0gMQp0bXBTTVszLCA4OjEzLCAyOjZdID0gMQp0bXBTTVs0LCA4OjEzLCAxNToxOV0gPSAxCnRtcFNNWzUsIDE1OjE5LCAyOjZdID0gMQp0bXBTTVs2LCAxNToxOSwgMTU6MTldID0gMQoKIyBQbG90IFNNcyBpbiBzaXggc3VicGxvdHMKcGFyKG1mcm93PWMoMiwgMykpCmZvciAoaSBpbiAxOjYpCnsKICBwbG90KHRtcFNNW2ksLF0sIG1haW49cGFzdGUoInRtcFNNIixpKSwgeGxhYj1OQSwgeWxhYj1OQSwgYm9yZGVyPU5BKQp9CmBgYAoKYGBge3IsIGVjaG89VFJVRSwgZmlnLmhlaWdodD0yLjUsIGZpZy53aWR0aD0yLjV9CiMgcGxvdCA2IHZlY3RvcmVkIFNNcyB0byBzaG93IGluZGVwZW5kZW5jZQpTTSA8LSBtYXRyaXgoYyhjKHRtcFNNWzEsLF0pLCBjKHRtcFNNWzIsLF0pLCBjKHRtcFNNWzMsLF0pLCBjKHRtcFNNWzQsLF0pLCBjKHRtcFNNWzUsLF0pLCBjKHRtcFNNWzYsLF0pKSwgNiwgMjEqMjEsIGJ5cm93PVRSVUUpCgpjb3JycGxvdChjb3IodChTTSkpLCBtYWluPSJcblNNIENvcnJlbGF0aW9uIEhlYXRtYXAiLCBtZXRob2Q9ImNvbG9yIikKYGBgCiMgUXVlc3Rpb24gMS40CmBgYHtyLCBlY2hvPVRSVUUsIGZpZy5oZWlnaHQ9MywgZmlnLndpZHRoPTV9CiMgQ29uc3RydWN0IEdhbW1hIHQgYW5kIEdhbW1hIHMKc2V0LnNlZWQoMTApCkdhbW1hX3QgPSBtYXRyaXgocm5vcm0oMjQwKjYsIG1lYW49MCwgc2Q9MC4yNSoqMC41KSwgMjQwLCA2KQpHYW1tYV9zID0gbWF0cml4KHJub3JtKDYqNDQxLCBtZWFuPTAsIHNkPTAuMDE1KiowLjUpLCA2LCA0NDEpCgojIFBsb3QgY29ycmVsYXRpb24KcGFyKG1mcm93PWMoMSwgMikpCmNvcnJwbG90KGNvcihHYW1tYV90KSwgbWFpbj0iXG5cblxuXG5UZW1wb3JhbCBOb2lzZSIsIAogICAgICAgICB4bGFiPU5BLCB5bGFiPU5BLCBtZXRob2Q9ImNvbG9yIikKY29ycnBsb3QoY29yKHQoR2FtbWFfcykpLCBtYWluPSJcblxuXG5cblNwYXRpYWwgTm9pc2UiLCAKICAgICAgICAgeGxhYj1OQSwgeWxhYj1OQSwgbWV0aG9kPSJjb2xvciIpCmBgYAoKYGBge3IsIGVjaG89VFJVRSwgZmlnLmhlaWdodD0yLCBmaWcud2lkdGg9NH0KZGYgPSBkYXRhLmZyYW1lKGMoR2FtbWFfdCkpCiMgUGxvdCBoaXN0b2dyYW0gb2YgYm90aCBub2lzZSBzb3VyY2UKZ2FtbWFfdF9wbG90IDwtIGdncGxvdChkZiwgYWVzKGMoYy5HYW1tYV90LikpKSsKICBnZW9tX2hpc3RvZ3JhbShhZXMoeT0uLmRlbnNpdHkuLiksCiAgYmlud2lkdGg9MC4xNSwgY29sb3I9ImRhcmtibHVlIiwgZmlsbD0ibGlnaHRibHVlIikgKwogIGxhYnModGl0bGU9IkRpc3RyaWJ1dGlvbiBvZiBUZW1wb3JhbCBOb2lzZSIseD0iVmFsdWUiLCB5ID0gIkRlbnNpdHkiKSArCiAgc3RhdF9mdW5jdGlvbihmdW4gPSBkbm9ybSwgYXJncyA9IGxpc3QobWVhbiA9IDAsIHNkID0gKDAuMjUpKiowLjUpKSArCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gKDEuOTYpKigwLjI1KSoqMC41LCBjb2xvcj0icmVkIikgKyAKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAtKDEuOTYpKigwLjI1KSoqMC41LCBjb2xvcj0icmVkIikKCmRmID0gZGF0YS5mcmFtZShjKEdhbW1hX3MpKQpnYW1tYV9zX3Bsb3QgPC1nZ3Bsb3QoZGYsIGFlcyhjKGMuR2FtbWFfcy4pKSkrCiAgZ2VvbV9oaXN0b2dyYW0oYWVzKHk9Li5kZW5zaXR5Li4pLAogIGJpbndpZHRoPTAuMDUsIGNvbG9yPSJkYXJrYmx1ZSIsIGZpbGw9ImxpZ2h0Ymx1ZSIpICsKICBsYWJzKHRpdGxlPSJEaXN0cmlidXRpb24gb2YgU3BhdGlhbCBOb2lzZSIseD0iVmFsdWUiLCB5ID0gIkRlbnNpdHkiKSArCiAgc3RhdF9mdW5jdGlvbihmdW4gPSBkbm9ybSwgYXJncyA9IGxpc3QobWVhbiA9IDAsIHNkID0gKDAuMDE1KSoqMC41KSkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9ICgxLjk2KSooMC4wMTUpKiowLjUsIGNvbG9yPSJyZWQiKSArIAogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IC0oMS45NikqKDAuMDE1KSoqMC41LCBjb2xvcj0icmVkIikKCnBsb3RfZ3JpZChnYW1tYV90X3Bsb3QsIGdhbW1hX3NfcGxvdCkKYGBgCgpgYGB7ciwgZWNobz1UUlVFLCBmaWcuaGVpZ2h0PTMsIGZpZy53aWR0aD0yLjV9CkdhbW1hX3RzIDwtIEdhbW1hX3QlKiVHYW1tYV9zCiMgQ2hlY2sgY29ycmVsYXRpb24gb2YgdGhlIHByb2R1Y3Qgb2YgR2FtbWEgdCBhbmQgR2FtbWEgcwpjb3JyIDwtIGNvcihHYW1tYV90cykKCiMgUmFuZG9tbHkgc2FtcGxlIDEwIHZhcmlhYmxlcyB0byB2aXN1YWxpc2UKc2FtcGxlcyA8LSBzYW1wbGUgKGMoMTo0NDEpLCBzaXplPTEwLCByZXBsYWNlPUYpCmNvcnJwbG90KGNvcnJbc2FtcGxlcywgc2FtcGxlc10sIG1ldGhvZD0iY29sb3IiLCBtYWluPSJcbkNvcnJlbGF0aW9uIGJldHdlZW4gNDQxIHZhcmlhYmxlcyAoU2FtcGxlZCkiLCB0eXBlID0gImxvd2VyIiwgdGwuc3J0ID0gMzApCmBgYAoKIyBRdWVzdGlvbiAxLjUKYGBge3J9CiMgR2VuZXJhdGUgYSBzeW50aGV0aWMgZGF0YXNldCBYClRDIDwtIGFzLm1hdHJpeChUQywgNiwgMjQwKQpYIDwtIChUQyArIEdhbW1hX3QpICUqJSAoU00gKyBHYW1tYV9zKQoKIyBDaGVjayBjb21wYXRpYmlsaXR5CmMoZGltKFRDKSwgZGltKEdhbW1hX3MpKQpjKGRpbShHYW1tYV90KSwgZGltKFNNKSkKZGltKFgpCmBgYAoKYGBge3J9CiMgUGxvdCBhdCBsZWFzdCAxMDAgcmFuZG9tbHkgc2VsZWN0ZWQgdGltZS1zZXJpZXMgZnJvbSBYCnNhbXBsZXMgPC0gZGF0YS5mcmFtZShuPTE6MjQwLCBYWywgc2FtcGxlLmludCgyNDAsIDEwMCldKQpzYW1wbGVzIDwtIG1lbHQoc2FtcGxlcywgaWQudmFycyA9ICJuIikKZ2dwbG90KHNhbXBsZXMsIGFlcyh4PW4sIHk9dmFsdWUsIGNvbD12YXJpYWJsZSkpICsgZ2VvbV9saW5lKCkgKyAKICBsYWJzKHRpdGxlPSIxMDAgUmFtZG9ubHkgU2VsZWN0ZWQgVGltZS1zZXJpZXMiLCB4PSIiLCB5PSJWYWx1ZSIpICsKICB0aGVtZShwbG90Lm1hcmdpbiA9IG1hcmdpbigxLjUsLjUsMS41LC41LCAiY20iKSwKICAgICAgICBsZWdlbmQua2V5LmhlaWdodD0gdW5pdCgwLjEsICdjbScpLAogICAgICAgIGxlZ2VuZC5rZXkud2lkdGg9IHVuaXQoMC4xLCAnY20nKSkKYGBgCgpgYGB7ciwgZWNobz1UUlVFLCBmaWcuaGVpZ2h0PTEuNSwgZmlnLndpZHRoPTIuNn0KIyBwbG90IHZhcmlhbmNlIG9mIGFsbCA0NDEgdmFyaWFibGVzCnZhcmlhbmNlcyA8LSBjb2xWYXJzKFgpCnZhcmlhbmNlcyA8LSBkYXRhLmZyYW1lKG49MTo0NDEsIHZhcmlhbmNlcykKZ2dwbG90KHZhcmlhbmNlcywgYWVzKHg9biwgeT12YXJpYW5jZXMpKSArZ2VvbV9wb2ludChjb2xvcj0ic3RlZWxibHVlIikgKyAKICBsYWJzKHRpdGxlPSJWYXJpYW5jZSBvZiA0NDEgdmFyaWFibGVzIiwgeD0iVmFyaWFibGVzIiwgeT0iVmFyaWFuY2VzIikKYGBgCmBgYHtyfQojIFN0YW5kYXJkaXplIFgKWCA9IHNjYWxlKFgpCmBgYAoKCiMgUXVlc3Rpb24gMi4xCmBgYHtyfQpEIDwtIGFzLm1hdHJpeChUQykKCiMgRXN0aW1hdGUgQSB1c2luZyBsZWFzdCBzcXVhcmUgZXN0aW1hdGlvbgpBX2xzciA8LSBhYnMoc29sdmUodChEKSUqJUQpJSoldChEKSUqJVgpCgojIENhbGN1bGF0ZSBEX2xzcgpEX2xzciA8LSBYJSoldChBX2xzcikKYGBgCgpgYGB7ciwgZWNobz1UUlVFLCBmaWcuaGVpZ2h0PTQsIGZpZy53aWR0aD0zLjV9CnBhbCA8LSBjb2xvclJhbXBQYWxldHRlKGMoImJsdWUiLCAieWVsbG93IikpCmxheW91dChtYXRyaXgoYygxLCAyLCAzLCA0LCA1LCA2LCA3LCA4LCA5LCAxMCwgMTEsIDEyKSwgNiwgMiwgbmNvbCA9IDIpLCAKICAgICAgIHdpZHRocz1jKDIuNSw1KSkKcGFyKG1hciA9IGMoMiw1LDEuNSwxLjUpKQpmb3IgKGkgaW4gMTo2KSAKewogIEFfbHNyX3AgPC0gQV9sc3JbaSxdCiAgRF9sc3JfcCA8LSBEX2xzclssaV0KICBkaW0oQV9sc3JfcCkgPC0gYygyMSwgMjEpCiAgcGxvdChBX2xzcl9wLCBib3JkZXI9RiwgY29sID0gcGFsKDQwKSwgeGxhYj1OQSwgeWxhYj1OQSwgCiAgICAgICBtYWluPXBhc3RlKCJSZXRyaWV2ZWQgU3BhdGlhbCBNYXAiLGkpKQogIAogIHBsb3QoRF9sc3JfcCwgdHlwZSA9ICJsIiwgeGxhYj1OQSwgeWxhYj1OQSwgCiAgICAgICBtYWluPXBhc3RlKCJSZXRyaWV2ZWQgVGltZSBDb3Vyc2UiLGkpKQp9CmBgYAoKYGBge3IsIGVjaG89VFJVRSwgZmlnLmhlaWdodD0xLjY1LCBmaWcud2lkdGg9NH0KZGYgPC0gZGF0YS5mcmFtZShEX2xzclssM10sIFhbLDMwXSkKY29sbmFtZXMoZGYpIDwtIGMoIkRfbHNyXzMiLCAiWF8zMCIpCnBsb3RfMSA8LSBnZ3Bsb3QoZGYsIGFlcyh5PURfbHNyXzMsIHg9WF8zMCkpICsgZ2VvbV9wb2ludChjb2xvciA9ICJzdGVlbGJsdWUiKSArCiAgbGFicyh0aXRsZT1leHByZXNzaW9uKHBhc3RlKERbbHNyXSwgIiAgY29sdW1uIDMgdnMgWCBjb2x1bW4gMzAiKSksCiAgICAgICB4PWV4cHJlc3Npb24ocGFzdGUoRFtsc3JdLCAiIGNvbHVtbiAzIikpLCB5ID0gIlggY29sdW1uIDMwIikgKyBnZW9tX3Ntb290aCgpCgpkZiA8LSBkYXRhLmZyYW1lKERfbHNyWyw0XSwgWFssMzBdKQpjb2xuYW1lcyhkZikgPC0gYygiRF9sc3JfNCIsICJYXzMwIikKcGxvdF8yIDwtIGdncGxvdChkZiwgYWVzKHk9RF9sc3JfNCwgeD1YXzMwKSkgKyBnZW9tX3BvaW50KGNvbG9yID0gInN0ZWVsYmx1ZSIpICsKICBsYWJzKHRpdGxlPWV4cHJlc3Npb24ocGFzdGUoRFtsc3JdLCAiICBjb2x1bW4gNCB2cyBYIGNvbHVtbiAzMCIpKSwKICAgICAgIHg9ZXhwcmVzc2lvbihwYXN0ZShEW2xzcl0sICIgY29sdW1uIDQiKSksIHkgPSAiWCBjb2x1bW4gMzAiKSArIGdlb21fc21vb3RoKCkKCnN1cHByZXNzTWVzc2FnZXMocGxvdF9ncmlkKHBsb3RfMSwgcGxvdF8yKSkKYGBgCiMgUXVlc3Rpb24gMi4yCmBgYHtyfQpsYW1iZGEgPC0gMC41CkkgPC0gbWF0cml4KDAsIDYsIDYpCmRpYWcoSSkgPC0gMQojIENvbXB1dGUgc2NhbGFyIHZhbHVlCnBlbmFsdHlfdGVybSA8LSBsYW1iZGEgKiBWCiMgRXN0aW1hdGUgQV9yciBhbmQgRF9ycgpBX3JyIDwtIGFicyhzb2x2ZSh0KEQpJSolRCArIHBlbmFsdHlfdGVybSpJKSUqJXQoRCklKiVYKQpEX3JyIDwtIFglKiV0KEFfcnIpCgpDX3Rsc3IgPSAwCkNfdHJyID0gMApmb3IgKGkgaW4gMTo2KQp7CiAgQ190bHNyW2ldID0gY29yKFRDWyxpXSwgRF9sc3JbLGldKQogIENfdHJyW2ldID0gY29yKFRDWyxpXSwgRF9yclssaV0pCn0KCiMgQ29tcGFyZSB0aGUgc3VtIG9mIHRoZXNlIHR3byBjb3JyZWxhdGlvbiB2ZWN0b3JzCnByaW50KHBhc3RlKCJTdW0gb2YgQ190bHNyOiIsIHN1bShDX3Rsc3IpKSkKcHJpbnQocGFzdGUoIlN1bSBvZiBDX3RycjoiLCBzdW0oQ190cnIpKSkKYGBgCgpgYGB7ciwgZWNobz1UUlVFLCBmaWcuaGVpZ2h0PTEuMywgZmlnLndpZHRoPTR9CiMgU2V0IGxhbWJkYSB2YWx1ZSBhcyAxMDAwIGZvciBwbG90IHZpc3VhbGlzYXRpb24KbGFtYmRhIDwtIDEwMDAKIyBDb21wdXRlIHNjYWxhciB2YWx1ZQpwZW5hbHR5X3Rlcm0gPC0gbGFtYmRhICogVgojIEVzdGltYXRlIEFfcnIgYW5kIERfcnIKQV9ycl8xMDAwIDwtIGFicyhzb2x2ZSh0KEQpJSolRCArIHBlbmFsdHlfdGVybSpJKSUqJXQoRCklKiVYKQoKIyBwbG90IGFuZCBjb21wYXJlIHRoZSBzaHJpbmthZ2Ugb2YgdmFyaWFibGVzCnZhbHVlcyA8LSBkYXRhLmZyYW1lKG49MTo0NDEsIHZhbD1BX2xzclsxLF0pCkFfbHNyX3Bsb3QgPC0gZ2dwbG90KHZhbHVlcywgYWVzKHg9biwgeT12YWwpKSArCiAgZ2VvbV9wb2ludChjb2xvciA9ICJzdGVlbGJsdWUiKSArIAogIGxhYnModGl0bGU9ZXhwcmVzc2lvbihBW2xzcl0pLCB5PSJWYWx1ZXMiLCB4PSJJbmRleCIpCgp2YWx1ZXMgPC0gZGF0YS5mcmFtZShuPTE6NDQxLCB2YWw9QV9ycl8xMDAwWzEsXSkKQV9ycl9wbG90IDwtIGdncGxvdCh2YWx1ZXMsIGFlcyh4PW4sIHk9dmFsKSkgKwogIGdlb21fcG9pbnQoY29sb3IgPSAic3RlZWxibHVlIikgKyAKICBsYWJzKHRpdGxlPWV4cHJlc3Npb24ocGFzdGUoQVtycl0sICIgICAiLCBsYW1iZGEsICIgPSAxMDAwIikpLCB5PSJWYWx1ZXMiLAogICAgICAgeD0iSW5kZXgiKQoKc3VwcHJlc3NNZXNzYWdlcyhwbG90X2dyaWQoQV9sc3JfcGxvdCwgQV9ycl9wbG90KSkKYGBgCiMgUXVlc3Rpb24gMi4zCmBgYHtyfQojIFNlbGVjdCDPgSBiZXR3ZWVuIDAgYW5kIDEgd2l0aCBhbiBpbnRlcnZhbCBvZiAwLjA1CnJob192YWx1ZXMgPC0gc2VxKGZyb209MCwgdG89MSwgYnk9MC4wNSkKTVNFIDwtIHJlcCgwLCAyMSkKYGBgCgpgYGB7cn0KZm9yIChqIGluIDE6MjEpCnsKICByaG8gPC0gcmhvX3ZhbHVlc1tqXQogIHN0ZXAgPC0gMS8obm9ybShUQyAlKiUgdChUQykpICogMS4xKQogIHRociA8LSByaG8qTipzdGVwCiAgQW8gPC0gbWF0cml4KDAsIG5zcmNzLCAxKQogIEEgPC0gbWF0cml4KDAsIG5zcmNzLCAxKQogIEFsciA8LSBtYXRyaXgoMCwgbnNyY3MsIHgxKngyKQogIE1TRV90ZW1wIDwtIHJlcCgwLCAxMCkKICAKICAjIEdlbmVyYXRlIG5ldyBzdGFuZGFyZGl6ZWQgWCBpbiB0ZXJtcyBvZiBHYW1tYSB0IGFuZCBHYW1tYSBzCiAgZm9yIChhIGluIDE6MTApIHsKICAgIHNldC5zZWVkKGEpCiAgICBHYW1tYV90X3RlbXAgPSBtYXRyaXgocm5vcm0oMjQwKjYsIG1lYW49MCwgc2Q9MC4yNSoqMC41KSwgMjQwLCA2KQogICAgR2FtbWFfc190ZW1wID0gbWF0cml4KHJub3JtKDYqNDQxLCBtZWFuPTAsIHNkPTAuMDE1KiowLjUpLCA2LCA0NDEpCiAgICBYX3RlbXAgPC0gKFRDICsgR2FtbWFfdF90ZW1wKSAlKiUgKFNNICsgR2FtbWFfc190ZW1wKQogICAgWF90ZW1wIDwtIHNjYWxlKFhfdGVtcCkKICAgIAogICAgZm9yIChrIGluIDE6KHgxKngyKSkgewogICAgICBBIDwtIEFvK3N0ZXAqKHQoVEMpICUqJSAoWF90ZW1wWyxrXS0oVEMlKiVBbykpKQogICAgICBBIDwtICgxLygxK3RocikpICogKHNpZ24oQSkqcG1heChyZXBsaWNhdGUobnNyY3MsIDApLCBhYnMoQSktdGhyKSkKICAgICAgCiAgICAgIGZvciAoaSBpbiAxOjEwKSB7CiAgICAgICAgQW8gPC0gQQogICAgICAgIEEgPC0gQW8rc3RlcCAqICh0KFRDKSUqJShYX3RlbXBbLGtdLShUQyUqJUFvKSkpCiAgICAgICAgQSA8LSAoMS8oMSt0aHIpKSAqIChzaWduKEEpKnBtYXgocmVwbGljYXRlKG5zcmNzLCAwKSwgYWJzKEEpLXRocikpCiAgICAgICAgfQogICAgICBBbHJbLGtdIDwtIEEKICAgIH0KICAgIEFsciA8LSBBbHIKICAgIERsciA8LSBYX3RlbXAlKiV0KEFscikKICAgIE1TRV90ZW1wW2FdIDwtIG5vcm0oKFhfdGVtcCAtIERsciUqJUFsciksIHR5cGU9J0YnKV4yLyhOKlYpCiAgfQogIAogICMgQ29tcHV0ZSB0aGUgYXZlcmFnZSBNU0Ugb2YgdGhlIGN1cnJlbnQgz4EKICBNU0Vbal0gPC0gbWVhbihNU0VfdGVtcCkKfQpgYGAKCmBgYHtyLCBlY2hvPVRSVUUsIGZpZy5oZWlnaHQ9MiwgZmlnLndpZHRoPTR9CiMgcGxvdCBhdmVyYWdlIG9mIE1TRSBvdmVyIHRoZXNlIDEwIHJlYWxpemF0aW9ucyBhZ2FpbnN0IGVhY2ggdmFsdWUgb2Ygz4EKCmRmIDwtIGRhdGEuZnJhbWUocmhvX3ZhbHVlcywgTVNFKQpwbG90XzEgPC0gZ2dwbG90KGRmLCBhZXMoeD1yaG9fdmFsdWVzLCB5PU1TRSkpICsgZ2VvbV9wb2ludChjb2xvcj0ic3RlZWxibHVlIikgKwogIGxhYnModGl0bGU9IkF2ZXJhZ2UgTVNFIGZvciBlYWNoIM+BIiwgeD0iz4EgdmFsdWVzIiwgeT0iQXZlcmFnZSBNU0UiKQoKZGYgPC0gZGZbMTA6MjEsXQpwbG90XzIgPC0gZ2dwbG90KGRmLCBhZXMoeD1yaG9fdmFsdWVzLCB5PU1TRSkpICsgZ2VvbV9wb2ludChjb2xvcj0ic3RlZWxibHVlIikgKwogIGxhYnModGl0bGU9IkF2ZXJhZ2UgTVNFIGZvciBlYWNoIM+BIChab29tZWQgSW4pIiwgeD0iz4EgdmFsdWVzIiwgeT0iQXZlcmFnZSBNU0UiKQoKc3VwcHJlc3NNZXNzYWdlcyhwbG90X2dyaWQocGxvdF8xLCBwbG90XzIpKQpgYGAKIyBRdWVzdGlvbiAyLjQKYGBge3J9CiMgRXN0aW1hdGUgQV9sciBhbmQgRF9sciB1c2luZyB0aGUgc2VsZWN0ZWQgz4EgdmFsdWUKcmhvIDwtIDAuNgpzdGVwIDwtIDEvKG5vcm0oVEMgJSolIHQoVEMpKSAqIDEuMSkKdGhyIDwtIHJobypOKnN0ZXAKQW8gPC0gbWF0cml4KDAsIG5zcmNzLCAxKQpBIDwtIG1hdHJpeCgwLCBuc3JjcywgMSkKQV9sciA8LSBtYXRyaXgoMCwgbnNyY3MsIHgxKngyKQoKZm9yIChrIGluIDE6KHgxKngyKSkgewogIEEgPC0gQW8rc3RlcCoodChUQykgJSolIChYWyxrXS0oVEMlKiVBbykpKQogIEEgPC0gKDEvKDErdGhyKSkgKiAoc2lnbihBKSpwbWF4KHJlcGxpY2F0ZShuc3JjcywgMCksIGFicyhBKS10aHIpKQogIAogIGZvciAoaSBpbiAxOjEwKSB7CiAgICBBbyA8LSBBCiAgICBBIDwtIEFvK3N0ZXAgKiAodChUQyklKiUoWFssa10tKFRDJSolQW8pKSkKICAgIEEgPC0gKDEvKDErdGhyKSkgKiAoc2lnbihBKSpwbWF4KHJlcGxpY2F0ZShuc3JjcywgMCksIGFicyhBKS10aHIpKQogIH0KICBBX2xyWyxrXSA8LSBBCn0KQV9sciA8LSBhYnMoQV9scikKRF9sciA8LSBYJSoldChBX2xyKQpgYGAKCmBgYHtyfQojIGVzdGltYXRlIGZvdXIgY29ycmVsYXRpb24gdmVjdG9ycwpDX3RyciA8LSAwCkNfc3JyIDwtIDAKQ190bHIgPC0gMApDX3NsciA8LSAwCgpmb3IgKGkgaW4gMTo2KQp7CiAgQ190cnJbaV0gPSBjb3IoVENbLGldLCBEX3JyWyxpXSkKICBDX3RscltpXSA9IGNvcihUQ1ssaV0sIERfbHJbLGldKQogIENfc3JyW2ldID0gY29yKFNNW2ksXSwgQV9ycltpLF0pCiAgQ19zbHJbaV0gPSBjb3IoU01baSxdLCBBX2xyW2ksXSkgIAp9CnN1bShDX3RycikKc3VtKENfdGxyKQpzdW0oQ19zcnIpCnN1bShDX3NscikKYGBgCgpgYGB7ciwgZWNobz1UUlVFLCBmaWcuaGVpZ2h0PTUsIGZpZy53aWR0aD01LjV9CmxheW91dChtYXRyaXgoc2VxKGZyb209MSwgdG89MjQsIGJ5PTEpLCA2LCA0LCBuY29sID0gNCksIAogICAgICAgd2lkdGhzPWMoMyw1LCAzLCA1KSkKcGFyKG1hciA9IGMoMyw1LDEuNSwyKSkKZm9yIChpIGluIDE6NikgewogIEFfcnJfcCA8LSBBX3JyW2ksXQogIEFfbHJfcCA8LSBBX2xyW2ksXQogIGRpbShBX3JyX3ApIDwtIGMoMjEsIDIxKQogIGRpbShBX2xyX3ApIDwtIGMoMjEsIDIxKQogIHBsb3QoQV9ycl9wLCBib3JkZXI9RiwgY29sID0gcGFsKDQwKSwgeGxhYj1OQSwgeWxhYj1OQSwgCiAgICAgICBjZXgubWFpbj0xLjYsIGNleC5sYWI9MS41LCBjZXguYXhpcz0xLjIsIG1haW49cGFzdGUoIkFyciIsaSkpCiAgcGxvdChEX3JyWyxpXSwgeGxhYj1OQSwgeWxhYj1OQSwgbWFpbj1wYXN0ZSgiRHJyIixpKSwgCiAgICAgICBjZXgubWFpbj0xLjYsIGNleC5sYWI9MS41LCBjZXguYXhpcz0xLjIsIGNvbD0iYmx1ZSIsIHR5cGU9ImwiKQogIHBsb3QoQV9scl9wLCBib3JkZXI9RiwgY29sID0gcGFsKDQwKSwgeGxhYj1OQSwgeWxhYj1OQSwgCiAgICAgICBjZXgubWFpbj0xLjYsIGNleC5sYWI9MS41LCBjZXguYXhpcz0xLjIsIG1haW49cGFzdGUoIkFsciIsaSkpCiAgcGxvdChEX2xyWyxpXSwgeGxhYj1OQSwgeWxhYj1OQSwgbWFpbj1wYXN0ZSgiRGxyIixpKSwKICAgICAgIGNleC5tYWluPTEuNiwgY2V4LmxhYj0xLjUsIGNleC5heGlzPTEuMiwgY29sPSJibHVlIiwgdHlwZT0ibCIpCn0KYGBgCiMgUXVlc3Rpb24gMi41IHVzaW5nIFNWRApgYGB7ciwgZWNobz1UUlVFLCBmaWcuaGVpZ2h0PTUsIGZpZy53aWR0aD02fQoKcGNyIDwtIHN2ZChEKQpaIDwtIHBjciR1CnBhcihtZnJvdz1jKDYsIDIpLCBtYXIgPSBjKDIsMywxLjUsMikpCgojIFBsb3QgdGhlIHJlZ3Jlc3NvcnMgaW4gWiBhbmQgc291cmNlIFRDcyBzaWRlIGJ5IHNpZGUKZm9yIChpIGluIDE6NikgewogIHBsb3QoVENbLGldLCB0eXBlPSJsIiwgbWFpbj1wYXN0ZSgiVEMiLCBpKSwgY2V4Lm1haW49MS41LCBjZXgubGFiPTEuNSwgY2V4LmF4aXM9MS41KQogIHBsb3QoWlssaV0sIHR5cGU9ImwiLCBtYWluPXBhc3RlKCJaIiwgaSksIGNleC5tYWluPTEuNSwgY2V4LmxhYj0xLjUsIGNleC5heGlzPTEuNSkKfQpgYGAKYGBge3IsIGVjaG89VFJVRSwgZmlnLmhlaWdodD0xLjUsIGZpZy53aWR0aD0zfQpldiA8LSBwY3IkZApldiA8LSBkYXRhLmZyYW1lKGluZGV4PXNlcSgxLCA2KSwgZXYpCnBsb3QoZXYsIGNvbD0iYmx1ZSIsIG1haW49IkVpZ2VudmFsdWVzIG9mIGVhY2ggUEMiLCB4bGFiPSJQQ3MiLCB5bGFiPSJFaWdlbnZhbHVlcyIsCiAgICAgY2V4Lm1haW49MS40LCBjZXgubGFiPTEuNCwgY2V4LmF4aXM9MS40KQpgYGAKCmBgYHtyfQojIEVzdGltYXRlIEFfcGNyIGFuZCBEX3BjciB1c2luZyB0aGUgz4EgPSAwLjAwMQptYXRyaXggPC0gWgpyaG8gPC0gMC4wMDEKc3RlcCA8LSAxLyhub3JtKFogJSolIHQoWikpICogMS4xKQp0aHIgPC0gcmhvKk4qc3RlcApBbyA8LSBtYXRyaXgoMCwgbnNyY3MsIDEpCkEgPC0gbWF0cml4KDAsIG5zcmNzLCAxKQpBX3BjciA8LSBtYXRyaXgoMCwgbnNyY3MsIHgxKngyKQoKZm9yIChrIGluIDE6KHgxKngyKSkgewogIEEgPC0gQW8rc3RlcCoodChtYXRyaXgpICUqJSAoWFssa10tKG1hdHJpeCUqJUFvKSkpCiAgQSA8LSAoMS8oMSt0aHIpKSAqIChzaWduKEEpKnBtYXgocmVwbGljYXRlKG5zcmNzLCAwKSwgYWJzKEEpLXRocikpCiAgCiAgZm9yIChpIGluIDE6MTApIHsKICAgIEFvIDwtIEEKICAgIEEgPC0gQW8rc3RlcCAqICh0KG1hdHJpeCklKiUoWFssa10tKG1hdHJpeCUqJUFvKSkpCiAgICBBIDwtICgxLygxK3RocikpICogKHNpZ24oQSkqcG1heChyZXBsaWNhdGUobnNyY3MsIDApLCBhYnMoQSktdGhyKSkKICB9CiAgQV9wY3JbLGtdIDwtIEEKfQpBX3BjciA8LSBhYnMoQV9wY3IpCkRfcGNyIDwtIFglKiV0KEFfcGNyKQpgYGAKCmBgYHtyLCBlY2hvPVRSVUUsIGZpZy5oZWlnaHQ9NCwgZmlnLndpZHRoPTMuNX0KbGF5b3V0KG1hdHJpeChjKDEsIDIsIDMsIDQsIDUsIDYsIDcsIDgsIDksIDEwLCAxMSwgMTIpLCA2LCAyLCBuY29sID0gMiksIAogICAgICAgd2lkdGhzPWMoMC44LDIpKQpwYXIobWFyID0gYygyLDUsMS41LDEuNSkpCgojIFBsb3QgdGhlIHJlc3VsdHMgb2YgRFBDUiBhbmQgQVBDUiBzaWRlIGJ5IHNpZGUKZm9yIChpIGluIDE6NikgCnsKICBBX3Bjcl9wIDwtIEFfcGNyW2ksXQogIERfcGNyX3AgPC0gRF9wY3JbLGldCiAgZGltKEFfcGNyX3ApIDwtIGMoMjEsIDIxKQogIHBsb3QoQV9wY3JfcCwgYm9yZGVyPUYsIGNvbCA9IHBhbCg0MCksIHhsYWI9TkEsIHlsYWI9TkEsIAogICAgICAgbWFpbj1wYXN0ZSgiQSBwY3IiLGkpKQogIAogIHBsb3QoRF9wY3JfcCwgdHlwZSA9ICJsIiwgeGxhYj1OQSwgeWxhYj1OQSwgCiAgICAgICBtYWluPXBhc3RlKCJEIHBjciIsaSkpCn0KYGBgCgojIFF1ZXN0aW9uIDIuNSB1c2luZyBQUkNPTVAKYGBge3IsIGVjaG89VFJVRSwgZmlnLmhlaWdodD0xLjUsIGZpZy53aWR0aD00fQptIDwtIHByY29tcChELCBjZW50ZXIgPSBUUlVFLCBzY2FsZS49VFJVRSkKWiA8LSBtJHgKcGV2IDwtIGZ2aXpfZWlnKG0pICsgCiAgbGFicyh0aXRsZT0iUGVyY2VudCBvZiBleHBsYWluZWQgdmFyaWFuY2VzIG9mIFBDcyIsIHg9IlBDcyIpCmVpZ2VudmFsdWVzIDwtIGZ2aXpfZWlnKG0sIGNob2ljZT0iZWlnZW52YWx1ZSIsIGdlb209ImxpbmUiKSArIAogIGxhYnModGl0bGU9IkVpZ2VudmFsdWUgb2YgZWFjaCBQQyIsIHg9IlBDUyIpCnN1cHByZXNzTWVzc2FnZXMocGxvdF9ncmlkKHBldiwgZWlnZW52YWx1ZXMpKQpgYGAKCmBgYHtyLCBlY2hvPVRSVUUsIGZpZy5oZWlnaHQ9NywgZmlnLndpZHRoPTZ9CnBhcihtZnJvdz1jKDYsIDIpLCBtYXIgPSBjKDIsMywxLjUsMikpCgojIFBsb3QgdGhlIHJlZ3Jlc3NvcnMgaW4gWiBhbmQgc291cmNlIFRDcyBzaWRlIGJ5IHNpZGUKZm9yIChpIGluIDE6NikgewogIHBsb3QoVENbLGldLCB0eXBlPSJsIiwgbWFpbj1wYXN0ZSgiVEMiLCBpKSwgY2V4Lm1haW49MS41LCBjZXgubGFiPTEuNSwgY2V4LmF4aXM9MS41KQogIHBsb3QoWlssaV0sIHR5cGU9ImwiLCBtYWluPXBhc3RlKCJaIiwgaSksIGNleC5tYWluPTEuNSwgY2V4LmxhYj0xLjUsIGNleC5heGlzPTEuNSkKfQpgYGAKCmBgYHtyfQojIEVzdGltYXRlIEFfcGNyIGFuZCBEX3BjciB1c2luZyB0aGUgz4EgPSAwLjAwMQptYXRyaXggPC0gWgpyaG8gPC0gMC4wMDEKc3RlcCA8LSAxLyhub3JtKFogJSolIHQoWikpICogMS4xKQp0aHIgPC0gcmhvKk4qc3RlcApBbyA8LSBtYXRyaXgoMCwgbnNyY3MsIDEpCkEgPC0gbWF0cml4KDAsIG5zcmNzLCAxKQpBX3BjciA8LSBtYXRyaXgoMCwgbnNyY3MsIHgxKngyKQoKZm9yIChrIGluIDE6KHgxKngyKSkgewogIEEgPC0gQW8rc3RlcCoodChtYXRyaXgpICUqJSAoWFssa10tKG1hdHJpeCUqJUFvKSkpCiAgQSA8LSAoMS8oMSt0aHIpKSAqIChzaWduKEEpKnBtYXgocmVwbGljYXRlKG5zcmNzLCAwKSwgYWJzKEEpLXRocikpCiAgCiAgZm9yIChpIGluIDE6MTApIHsKICAgIEFvIDwtIEEKICAgIEEgPC0gQW8rc3RlcCAqICh0KG1hdHJpeCklKiUoWFssa10tKG1hdHJpeCUqJUFvKSkpCiAgICBBIDwtICgxLygxK3RocikpICogKHNpZ24oQSkqcG1heChyZXBsaWNhdGUobnNyY3MsIDApLCBhYnMoQSktdGhyKSkKICB9CiAgQV9wY3JbLGtdIDwtIEEKfQpBX3BjciA8LSBhYnMoQV9wY3IpCkRfcGNyIDwtIFglKiV0KEFfcGNyKQpgYGAKCmBgYHtyLCBlY2hvPVRSVUUsIGZpZy5oZWlnaHQ9NCwgZmlnLndpZHRoPTMuNX0KbGF5b3V0KG1hdHJpeChjKDEsIDIsIDMsIDQsIDUsIDYsIDcsIDgsIDksIDEwLCAxMSwgMTIpLCA2LCAyLCBuY29sID0gMiksIAogICAgICAgd2lkdGhzPWMoMC44LDIpKQpwYXIobWFyID0gYygyLDUsMS41LDEuNSkpCgojIFBsb3QgdGhlIHJlc3VsdHMgb2YgRFBDUiBhbmQgQVBDUiBzaWRlIGJ5IHNpZGUKZm9yIChpIGluIDE6NikgCnsKICBBX3Bjcl9wIDwtIEFfcGNyW2ksXQogIERfcGNyX3AgPC0gRF9wY3JbLGldCiAgZGltKEFfcGNyX3ApIDwtIGMoMjEsIDIxKQogIHBsb3QoQV9wY3JfcCwgYm9yZGVyPUYsIGNvbCA9IHBhbCg0MCksIHhsYWI9TkEsIHlsYWI9TkEsIAogICAgICAgbWFpbj1wYXN0ZSgiQSBwY3IiLGkpKQogIAogIHBsb3QoRF9wY3JfcCwgdHlwZSA9ICJsIiwgeGxhYj1OQSwgeWxhYj1OQSwgCiAgICAgICBtYWluPXBhc3RlKCJEIHBjciIsaSkpCn0KYGBgCgpgYGB7cn0KZGYgPC0gZGF0YS5mcmFtZShtJHJvdGF0aW9uKQpkZgpgYGAKCg==